Diffusion of hydrogen in graphite: A molecular dynamics simulation 
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Diffusion of atomic and molecular hydrogen in the interstitial space between graphite sheets 
has been studied by molecular dynamics simulations. Interatomic interactions were modeled by a 
tight-binding potential fitted to density-functional calculations. Atomic hydrogen is found to be 
bounded to C atoms, and its diffusion consists in jumping from a C atom to a neighboring one, with 
an activation energy of about 0.4 eV. Molecular hydrogen is less attached to the host sheets and 
diffuses faster than isolated H. At temperatures lower than 500 K, H2 diffuses with an activation 
energy of 89 meV, whereas at higher T its diffusion is enhanced by longer jumps of the molecule as 
well as by correlations between successive hops, yielding an effective activation energy of 190 meV. 

PACS numbers: 68.43.Jk, 68.43. Fg, 81.05.Uf 
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I. INTRODUCTION 

The past few years have seen extraordinary progress 
in the knowledge of carbon-based materials, and in par- 
ticular on those formed by C atoms displaying sp 2 hy- 
bridization. This is the case of materials discovered in 
last decades, such as carbon nanotubes and fullerenes, or 
in last years, such as graphene [HIH, apart from the more 
traditional graphite. 

Carbon-based systems, in general, are considered 
as candidates for hydrogen storage [1, 0]. More- 
over, chemisorption on two-dimensional systems, such as 
graphene or graphite surfaces, can be important for cat- 
alytic processes [5| . The interest on hydrogen as an impu- 
rity in solids and on surfaces is not new, and dates back 
to several decades. This is in principle one of the sim- 
plest impurities, but a deep understanding of its physical 
properties is complex due to its low mass, and requires 
the combination of advanced experimental and theoreti- 
cal methods [1, @] • In addition to its basic interest as an 
impurity, a relevant characteristic of hydrogen in solids 
and surfaces is its ability to form complexes and passi- 
vate defects, which has been extensively studied in the 
last thirty years 

Experimental investigations on atomic, isolated hydro- 
gen in graphite have been so far scarce, due to the dif- 
ficulty in detecting this impurity. Moreover, this prob- 
lem is complicated by the presence of a large amount 
of hydrogen trapped at the boundaries of graphite crys- 
tallites |9hT1|]. The stable hydrogen configurations in 
the bulk of graphite have been investigated in various 
theoretical works [T2l4l4j |. with particular emphasis on 
atomic and molecular forms of this impurity. Also, the 
diffusion, trapping, and recombination of hydrogen on a 
graphite surface have been studied by theoretical tech- 
niques |15l - ll8l |. In connection with this, atomic hydrogen 
on graphene has been investigated by several authors us- 
ing ab-initio methods [El [l9l - l23l | . It is generally accepted 
that chemisorption of a single hydrogen atom leads to the 
appearance of a defect-induced magnetic moment on the 
graphene sheet, along with a large structural distortion 



[19h21| ]. Recently, Ohldag et al. [24[ have discussed the 
role of hydrogen in room-temperature ferromagnetism at 
graphite surfaces from an x-ray dichroism analysis. 

For the storage of hydrogen in graphite one should 
consider, apart from atomic hydrogen, the presence 
of H2 molecules between the graphite sheets, which 
are ex pec ted to be physisorbed in the interlayer space 
d [H 13 • To study the diffusion of H 2 in the bulk of 
graphite, one has to take into account that most of the 
experimental measurements detect in fact the molecular 
diffusion through crystallite boundaries so that it is 
hard to obtain direct insight into the molecular diffusion 
in the interlayer space. 

In this paper, molecular dynamics (MD) simulations 
are used to investigate the diffusion of atomic and molec- 
ular hydrogen in graphite. Recently, the diffusion of 
atomic hydrogen on an isolated graphene sheet has been 
studied by a combination of path-integral molecular dy- 
namics simulations and transition-state theory, with spe- 
cial emphasis upon the appearance of quantum effects 
[2f|. Here, we will deal with classical molecular dynam- 
ics simulations, at temperatures high enough that such 
quantum effects become unimportant. This allows us 
to calculate directly diffusion coefficients, as the time in 
the classical simulations corresponds to real time, con- 
trary to path-integral simulations, where the simulation 
time does not strictly correspond to actual time, but ap- 
pears as a convenient computational way to derive time- 
independent thermodynamic properties. In the simula- 
tions presented here, the interatomic interactions have 
been modeled by a tight-binding (TB) potential fitted 
to density-functional calculations. The thermal behav- 
ior of hydrogen between graphite layers, as well as its 
diffusion in porous graphite have been addressed ear- 
lier by using MD simulations QiJ [3. This computa- 
tional technique has been also applied to study several 
pro pert ies of hydrogen in various carbon-based materi- 
als [26V 29] . In general, finite-temperature properties of 
hydrogen-related defects in solids have been investigated 
by ab-initio and TB molecular dynamics simulations (30l — 
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The paper is organized as follows. In Sec. II, we de- 
scribe the computational method employed in our cal- 
culations. Our results are presented in Sec. Ill, dealing 
with the diffusion of atomic and molecular hydrogen. In 
Sec. IV we summarize the main results. 



II. COMPUTATIONAL METHOD 

An important issue of the MD method is the adequate 
description of interatomic interactions, which should be 
as realistic as possible. Since effective classical potentials 
present many limitations to reproduce the many-body 
energy surface, specially in those situations where inter- 
atomic bonds may be either broken or formed, one should 
resort to self-consistent quantum-mechanical methods. 
However, density functional or Hartree-Fock based self- 
consistent potentials require computer resources that 
would restrict enormously the size of our simulation cell 
and/or the number of simulation steps. We found a 
compromise by employing an efficient tight-binding ef- 
fective Hamiltonian, based on density functional (DF) 
calculations [33j . The ability of TB methods to repro- 
duce different properties of solids and molecules has been 
reviewed by Goringe et al. [34| We checked the capa- 
bility of this DF-TB potential to predict frequencies of 
C-H vibrations. In particular, for a methane molecule 
it predicts in a harmonic approximation frequencies of 
3100 and 3242 cm -1 for C-H modes with symmetry A\ 
and T2, respectively [35[, to be compared with experi- 
mental values of 2917 and 3019 cm -1 I3Q|. Consider- 



ing the usual anharmonic shift (towards lower frequen- 
cies) associated to these modes, the agreement is sat- 
isfactory. A detailed analysis of vibrational frequencies 
of hydrocarbon molecules derived with the present DF- 
TB potential, including anharmonicities, can be found 
elsewhere [37], HI]. We have employed earlier this TB 
Hamiltonian to describe hydrogen-carbon interactions in 
diamond [HI, [3j| and graphene [25| . The TB energy con- 
sists of two parts, the first one is the sum of energies 
of occupied one-electron states, and the second one is 
given by a pairwise repulsive interatomic potential [33| . 
Since a correct description of the hydrogen molecule is 
essential for our purposes, special care was taken with 
the H-H pair potential, which has been taken as in our 
earlier study of molecular hydrogen in the silicon bulk 
[icj . This pair potential reproduces the main features of 
known effective interatomic potentials for H2 , such as the 
Morse potential fill ]. 

Molecular dynamics simulations were carried out in 
the NVE ensemble on a graphite supercell containing 
64 C atoms and one impurity (H or H2), and periodic 
boundary conditions were assumed. The simulation cell 
includes two graphite sheets, each one being a 4 x 4 
graphene supercell of size 4a = 9.84 A. We considered 
an AB layer stacking, so that both sheets are disposed in 
such a way that the center of each hexagonal ring of one 
of them lies over a C atom of the adjacent sheet. To hold 



this kind of stacking along a simulation run, thus avoid- 
ing diffusion of the graphite layers, the center-of-gravity 
of each layer was not allowed to move on the (x, y) plane. 
Note that in the following we will refer to the z direction 
as the one perpendicular to the graphite layers. The av- 
erage distance between sheets is a half of the supercell 
parameter along the z axis, and was taken to be 3.35 A. 
For the reciprocal-space sampling we have employed only 
the T point (k = 0), as the main effect of using a larger 
k set is a nearly constant shift in the total energy, with 
little influence in the calculation of energy differences. 

Sampling of the configuration space has been per- 
formed at temperatures between 300 and 2000 K. For 
a given temperature, a typical run consisted of 3 x 10 4 
MD steps for system equilibration, followed by 4 x 10 6 
steps for the calculation of ensemble average properties. 
In some cases, specially at low temperatures, we carried 
out longer simulations to reduce the statistical errors. In 
particular, for H2 at T < 400 K the simulations included 
8 x 10 6 MD steps. The equations of motion were inte- 
grated by using the standard Verlet algorithm [42[ . The 
time step At was taken in the range between 0.2 and 
0.5 fs, which was found to be appropriate for the atomic 
masses and temperatures studied here. Some checks were 
carried out for smaller values of At, yielding within er- 
ror bars the same diffusion coefficients as those reported 
below. The diffusion coefficient at different temperatures 
has been calculated from the long-time behavior of the 
mean-square displacement of the mobile species (H or 
H2) in the interlayer space. Although motion along the 
z coordinate (perpendicular to the graphite sheets) is al- 
lowed, it does not contribute to the long-time displace- 
ments, which are basically two-dimensional. Thus, the 
diffusion coefficient is given by: 



D 



1 {Ax{t)f + {Ay{t)f 
— am 

4 t^oo t 



(1) 



where the mean-square displacement in the x coordinate 
is 



(Ax(t)) 2 = ((x(t + t ) ~x(t )f 



(2) 



and a similar expression holds for (Ay(t)) 2 . Here (...) 
means an average over different zero times to along a 
MD trajectory. For molecular hydrogen, the coordinates 
x and y correspond to the center-of-mass of the molecule. 



III. RESULTS 

A. Atomic hydrogen 

We first discuss the lowest-energy configuration for 
atomic hydrogen in graphite, as derived from calcula- 
tions at T = K. The impurity binds to a C atom, which 
relaxes out of the sheet plane by 0.34 A, with a bond dis- 
tance between C and H of 1.17 A. This result is in line 
with those reported in the literature, with the breaking 
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of a 7r bond and the creation of an additional a bond, 
changing the hybridization of the involved C atom from 
sp 2 to sp 3 0, El d M, [2l|. This configuration with H 
attached to a C atom, which strongly relaxes off the sheet 
plane, is similar to that found for hydrogen adsorbed on 
an isolated graphene sheet SUM 

For its relation with the diffusion process, we have cal- 
culated the energy barrier to break a C-H bond and at- 
tach the hydrogen to a nearest C atom in the same sheet, 
and found an energy of 0.40 eV. Note that this value 
includes contributions of the appreciable relaxations of 
both carbon atoms involved in the whole process. It is 
interesting that we find the same energy barrier (within 
the precision of our method) for hydrogen to detach from 
a C atom in one sheet and jump to an opposite C atom 
in the adjacent sheet. This is not strange if one takes 
into account that a hydrogen atom linked to a C atom in 
a graphite sheet is at a distance of around 1.51 A from 
the sheet plane, a value close to half the distance be- 
tween graphite layers (1.67 A). Since H moves towards 
the middle plane to brake a C-H bond, the barriers for 
jumping to a C atom in the same or in the adjacent sheet 
are similar. The energy value found here is close to those 
obtained in earlier DF calculations: a value of 0.48 eV 
was reported in Ref. [l5[ , and an estimated barrier of 0.4- 
0.5 eV in Ref. [l3j]. We note that this barrier is clearly 
lower than the one found for hydrogen jumps on an iso- 
lated graphene layer, for which the interaction potential 
employed here gives an energy of 0.78 eV. 

We now turn to the results of our MD simulations at 
finite temperatures. From the calculations at T = K, 
one can expect that hydrogen will diffuse between the 
graphite sheets in a step-like fashion, dissociating from 
a C atom and attaching to a nearby one, in the same 
or in the adjacent sheet. This will cause an almost two- 
dimensional diffusion in the interlayer space. The likeli- 
hood of hydrogen jumping from an interlayer region to an 
adjacent one, crossing a graphite sheet is very low, since 
H has to climb an energy barrier of about 5 eV through 
an hexagonal ring of C atoms. 

In Fig. 1 we show the evolution of the z coordinate of 
hydrogen (direction perpendicular to the graphite sheets) 
along a MD simulation consisting of 10 6 steps (amount- 
ing to a time of 300 ps), at a temperature of 600 K. In 
spite of the large fluctuations in the coordinate z, one ob- 
serves that the hydrogen atom is mainly located on one 
of two planes at distance of about 1.5 or 1.8 A from the 
lower graphite layer. These two distances correspond to 
H linked to C atoms in the lower and upper sheets, re- 
spectively. In fact, from the calculations at zero temper- 
ature, we found that the stable sites for hydrogen (bound 
to C atoms) lie at 1.51 A from a graphite sheet, and at 
3.35 A - 1.51 A = 1.84 A from the adjacent one (see 
dashed lines in Fig. 1). In general, we observe along 
the MD simulations that the hydrogen jumps are rather 
short and direct, breaking a C-H bond and forming an- 
other one with a nearby carbon atom in the same or in 
an adjacent sheet. 



-< 



C3 



o 
o 
o 




100 200 

Simulation time (ps) 



FIG. 1: Coordinate z of hydrogen along a simulation run at 
600 K. The data shown include 10 6 MD steps, corresponding 
to a simulation time of 300 ps. The coordinate z around 1.5 
and 1.8 A corresponds to hydrogen attached to C atoms either 
in the lower or in the upper graphite sheet, z is measured from 
the lower sheet. 
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FIG. 2: Diffusion coefficient of atomic hydrogen in graphite, 
shown in an Arrhenius plot vs the inverse temperature. Error 
bars are on the order of the symbol size. The dashed line 
is a least-square fit to the data points, giving an effective 
activation energy E a — 0.38 eV. 



From the mean-square displacement of hydrogen along 
the MD simulations, we have calculated the diffusion co- 
efficient D in the interlayer space by using Eq. ([T]). The 
results are presented in Fig. 2 as a function of the inverse 
temperature in an Arrhenius plot. They can be fitted 
well to an expression D = Doexp(—E a /ltBT), with an 
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activation energy E a — 0.38 eV and a preexponential 
factor D = 2.0 x 10~ 3 cm 2 /s. This value corresponds 
to the order of magnitude expected for H diffusion in 
graphite. In fact, Dq can be expressed by the simple 
expression: Dq = (Ar) 2 z/o/4, where Ar is the distance 
between nearest adsorption sites (C atoms), and vq is a 
typical "attempt frequency" (the factor of four in the de- 
nominator takes into account the fact that the diffusion is 
two-dimensional, as in Eq. Taking Ar = 1.4 A and 
an attempt frequency corresponding to a C-H stretching 
mode (about 3000 cm -1 ), we find an estimation for the 
preexponential factor Dq — 4.4 x 10~ 3 cm 2 /s. This esti- 
mation is of the order of Dq derived from the simulations, 
but its numerical value is somewhat larger. This is not 
strange if one considers that hops of hydrogen from one 
layer to the adjacent one can contribute Ar = 1.4 A in the 
[x, y) plane, but also Ar = (jumping between opposite 
C atoms in the AB layer stacking), without contributing 
to the overall diffusion. 

At 1000 K we obtain from the MD simulations a rel- 
atively high diffusion coefficient D = 2.4 x 10~ 5 cm 2 /s. 
In practice, we can determine rather accurately the value 
of D down to temperatures in the order of 600 K, where 
D ~ 10~ 6 cm 2 /s. At lower T, the hydrogen diffusion 
is too slow to allow for a precise determination of D 
from the mean-square displacements. In fact at room 
temperature (T = 300 K) hydrogen jumps are observed 
very rarely in the simulations. In any case, an extrapo- 
lation of the fit shown in Fig. 2 yields at 300 K a value 
D = 9.2 x 10~ 10 cm 2 /s. The activation energy derived 
from the slope of the dashed line displayed in Fig. 2 (E a 
= 0.38 eV) is very close to the energy barrier calculated 
for breaking a C-H bond and linking the hydrogen to 
another C atom. This confirms the diffusion mechanism 
of atomic hydrogen expected from the calculations at T 
= K. 

At T lower than room temperature one expects the 
appearance of appreciable quantum effects in the hydro- 
gen diffusion. Such quantum effects cause an effective 
renormalization of the diffusion barrier, which is reduced 
with respect to the high-temperature value. This can be 
studied by a combination of transition-state theory with 
quantum path-integral simulations, as has been done for 
hydrogen on agraphene sheet [25| and in bulk semicon- 
ductors [1^, l43j | but is out of the scope of the present 
work. 



B. Molecular hydrogen 

Molecular hydrogen is expected to diffuse in graphite 
faster than atomic hydrogen, since the former will be 
less bonded to the host sheets, and therefore more mo- 
bile in the interlayer space. For an H2 molecule we find a 
minimum-energy position at an interstitial site between 
a carbon atom in a graphite sheet and an hexagonal ring 
in an adjacent sheet. At this position, the preferred 
orientation of the molecule is parallel to the graphite 
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FIG. 3: Probability distribution f(cp) for the angle cp between 
the H-H direction and the graphite sheets, as derived from 
MD simulations at three different temperatures: 300 K (solid 
line), 800 K (dashed line), and 1300 K (dotted line). 



planes, in agreement with earlier calculations based on 
density-functional theory (l3| . At finite temperatures 
the molecule will explore other positions and orientations 
with respect to the graphite layers. We define the angular 
probability distribution f(cp) of the angle cp between the 
H-H direction and the sheet plane such that the proba- 
bility P(tpi,if2) of observing an angle cp in the interval 
(^1,^2) is given by 



P(<Pl,<P2) 



f(cp) cosepdep 



(3) 



•Pi 



(i.e., cos cp takes into account the degeneracy of angle cp). 
In Fig. 3 we present the probability distribution f(cp), 
as derived from our molecular dynamics simulations at 
three different temperatures: 300 K (solid line), 800 K 
(dashed line), and 1300 K(dotted line). The distribution 
has a maximum at cp = (H-H parallel to the layers), 
and reaches its minimum for H-H perpendicular to the 
sheet plane (cp = 90°). As temperature increases, the 
probability distribution broadens slightly, but it remains 
as a peak centered at cp — 0. 

From the mean-square displacement of H2 along the 
molecular dynamics simulations we have calculated its 
diffusion coefficient at several temperatures. The results 
are presented in Fig. 4 as a function of the inverse tem- 
perature. One observes first that the diffusion coefficient 
of H2 is larger than that of atomic hydrogen in graphite, 
shown in Fig. 2. Thus, at 1000 K we find D = 2.5 x 10~ 4 
cm 2 /s for H2 to be compared with 2.4 x 10 -5 cm 2 /s for 
H. This difference is much larger at 300 K, since we find 
D — 6.9 x 10 -6 cm 2 /s for H2 vs the extrapolated value of 
D — 9.2 x 10~ 10 cm 2 /s for atomic hydrogen (about four 
orders of magnitude less) . 
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FIG. 4: Diffusion coefficient of molecular hydrogen in 
graphite as a function of the inverse temperature, as derived 
from our MD simulations. Squares represent data obtained 
from simulations where the C atoms of graphite were allowed 
to move, whereas circles correspond to simulations where the 
C atoms were kept fixed on their ideal equilibrium sites in 
the graphite sheets. Data shown as squares were fitted in- 
dependently at high and low temperatures, giving activation 
energies of 190 and 98 meV for T > 500 K and T < 500 K, 
respectively. Error bars of the simulation results are on the 
order of the symbol size. 



We note that the results presented for the diffusion co- 
efficient of H2 were derived from MD simulations where 
the graphite sheets are flexible, i.e. the C atoms are 
free to move along the simulations, and in particular the 
sheets can relax and adapt to the presence of the hydro- 
gen molecule. However, due to the lack of direct bonding 
between molecular hydrogen and graphite, one can ask it 
relaxation of the sheets does in fact affect the molecular 
diffusion. To obtain insight into this question we have 
carried out some MD simulations in which the C atoms 
were kept fixed at their unrelaxed (ideal) positions, and 
calculated the diffusion coefficient of H2. The results of 
these calculations are displayed in Fig. 4 as circles at 
temperatures higher than 1000 K. One notices that, at a 
given T, the value of D derived in this way is clearly lower 
than that found when the C atoms are allowed to move 
and relax in the presence of H2. In fact, at 1000 K the 
diffusion coefficient with fixed C atoms is smaller than 
10~ 6 cm 2 /s, the lower limit for a reliable determination 
of D in our simulations. All this indicates that motion 
and relaxation of the graphite sheets directly affects the 
diffusion of molecular hydrogen in graphite. 

A remarkable aspect of Fig. 4 is that the diffusion co- 
efficient of molecular hydrogen does not follow a single 
straight line in the Arrhenius plot, i.e. it does not dis- 
play a dependence of the form D = Dq exp(— Ea/ksT) in 
the whole temperature range considered here. It seems 
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FIG. 5: Mean residence time of molecular hydrogen in an 
interlayer site of graphite plotted vs the inverse temperature. 
Error bars are in the order of the symbol size. A least-square 
fit to the data points (dashed line) gives an activation energy 
of 97 meV. 



rather that one can distinguish two temperature regions 
with different effective activation energy E a . In fact, at 
T > 500 K those data can be fitted well with E a = 0.19 
eV, whereas at lower temperature E a seems to be about 
a factor of two lower (E a = 0.098 cV). 

Given the change in slope of the results presented in 
the Arrhenius plot of Fig. 4 (squares), we wonder whether 
at high and low temperatures the diffusion mechanism is 
different, or maybe the molecular jump rate from site to 
site in the interlayer space suffers some change apart from 
that predicted by a single activation energy. To clarify 
this question we have calculated the mean residence time 
of H2 on the interstitial sites along the MD simulations. 
To this end we have counted the number of molecular 
jumps in the molecular dynamics trajectories at different 
temperatures. The results for the mean residence time 
t are plotted in Fig. 5 vs the inverse temperature, t 
is found to follow a dependence t = tq exp(_E a //csT) in 
the whole temperature range under consideration, with 
an activation energy E a = 97 meV. This activation en- 
ergy coincides, within statistical error, with that derived 
from the diffusion coefficient of molecular hydrogen at 
temperatures lower than 500 K (E a = 98 meV). 

The available sites for H2 in the interlayer space of 
graphite form a honeycomb lattice, in which each site has 
three nearest neighbors. At low temperature (T < 500 
K) the molecular diffusion proceeds via jumps between 
nearest-neighbor adsorption sites in a random walk of 
the H 2 molecule in the interlayer space. At high temper- 
atures, however, we observe two mechanisms that con- 
tribute to enhance the molecular diffusion. First, we 
detected jumps longer than the distance between near- 
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est interstitial adsorption sites. In particular, we found 
direct jumps to second (next-nearest) and third neigh- 
boring interlayer sites, at a distance of 1.7 and 2 times 
that between nearest neighbors. At 1000 K, the fraction 
of these longer jumps amounts to about 9% of the total 
number of molecular jumps, versus around 1% at 500 K, 
and a negligible quantity of 300 K. 

At T > 500 K we also observed that the molecular 
jumps between nearest adsorption sites are correlated, in 
the sense that after a given jump, the probability that 
the next one occurs in the forward direction (the same 
direction as the previous one) is larger than in any other 
direction. This correlation contributes to enhance the 
diffusion of H2 , and at high temperatures the molecular 
diffusion cannot be considered as a random walk over the 
available adsorption sites. These two mechanisms (longer 
and correlated jumps) cause an increase in the diffusion 
coefficient, as compared with the behavior expected from 
an extrapolation of the low-temperature results. Such an 
increase translates into an enhancement of the slope in 
the Arrhenius plot of the diffusion coefficient, and conse- 
quently in the observation of an apparently larger effec- 
tive activation energy E a . This behavior illustrates the 
limitations of the interpretation of diffusion coefficients 
in Arrhenius plots. In fact, any temperature-dependent 
effect in the preexponential factor cannot be directly dis- 
tinguished in an Arrhenius plot from changes in the acti- 
vation energy. Our results of the mean residence time of 
H 2 in graphite provide evidence that the energy barrier 
for molecular motion remains constant in the tempera- 
ture range from 250 to 1600 K. 

IV. SUMMARY 

The main advantage of this kind of MD simulations of 
hydrogen in graphite is the possibility of calculating dif- 
fusion coefficients in a large range of temperatures, using 
an interatomic potential fitted to ab-initio calculations. 
Due to the large relaxation of the nearest C atoms, the 
migration of atomic hydrogen in graphite requires impor- 
tant motion of these atoms. Then, a hydrogen jump has 



to be viewed as a cooperative process involving a coupled 
motion of the impurity and the nearest host atoms. This 
effect is not a requirement for the diffusion of H2 , but in 
practice relaxation of the C atoms in the nearest graphite 
layers helps to enhance appreciably the molecular diffu- 
sion. This has been shown in Fig. 4 by comparing MD 
simulations with fixed or mobile carbon atoms. In fact, 
at 1000 K the diffusion coefficient is more than 100 times 
larger when the C atoms are allowed to relax out of their 
ideal positions. 

At 1000 K we obtained for molecular hydrogen a dif- 
fusion coefficient one order of magnitude larger than for 
atomic hydrogen. In fact, we found D = 2.5 x 10~ 4 
cm 2 /s for H 2 versus 2.4 x 10~ 5 cm 2 /s for H. This differ- 
ence increases to about four orders of magnitude at 300 
K: D = 6.9 x 10~ 6 cm 2 /s for H 2 , to be compared with 
an extrapolated value of D = 9.2 x 10~ 10 cm 2 /s 

The diffusion coefficients derived here for atomic and 
molecular hydrogen in graphite have to be considered as 
higher limits for the actually measured values. In fact, 
point or extended defects in the host layers will act as 
hydrogen traps, causing an appreciable reduction in the 
diffusion. More important, for a comparison with ex- 
perimental measurements, one has to take into account 
that an important part of the hydrogen is concentrated 
on the crystallite boundaries and in the micro-voids be- 
tween graphite granules. 

An interesting question is the dissociation of molecular 
hydrogen and its recombination in the interlayer space of 
graphite. A study of the combination of these processes 
with atomic and molecular diffusion lies, however, out of 
the scope of the present paper and remains as a challenge 
for future research. 
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